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ABSTRACT 


A  model  of  a  distributed  feedback  (DFB)  laser  was  implemented  in  MATLAB  and 
SIMULINK.  Using  the  laser  rate  equation,  the  model  was  simulated  to  obtain  general 
characteristics  of  the  chirp  of  the  laser’s  frequency.  The  simulations  were  controlled  by 
using  different  drive  current  waveforms,  based  on  various  bit-patterns,  data  rates,  and 
drive  current  values  (threshold  current  and  the  extinction  ratio).  Once  created,  the  laser 
drive  current  was  passed  to  the  SIMULINK  DFB  laser  model.  The  output  of  a  simulation 
provided  frequency  chirp,  laser  power  emitted,  photon  density,  and  carrier  density  data 

Two  sets  of  simulations  were  conducted.  The  first  set  of  simulations  focused  on 
the  data  rates  and  bit-patterns.  From  these  simulations  it  was  determined  that  the  transition 
from  a  ZERO  bit  to  a  ONE  bit  caused  the  greatest  frequency  excursions.  Also,  as  the  data 
rate  increases  the  maximum  frequency  excursion  increases.  Finally,  the  first  set  of 
simulations  revealed  that  the  predictability  of  the  chirp  decreases  as  the  data  rate  increases 

and  as  the  complexity  of  the  bit  pattern  increases. 

The  second  set  of  simulations  examined  the  effect  of  the  extinction  ratio  on 
frequency  chirp.  By  plotting  the  maximum  frequency  excursion  against  its  respective 
extinction  ratio,  it  was  determined  that  in  some  cases  the  maximum  frequency  excursions 
in  a  system  could  be  minimized. 
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I.  INTRODUCTION 


Wavelength-division  multiplexed  (WDM)  or  wavelength-divison  multiple  access 
(WDMA)  systems  have  been  proposed  to  increase  the  aggregate  bit  rates  of  fiber  optic 
systems.  Single  channel  bit  rates  as  high  as  10  Gb/s  are  envisioned  to  support  networks 
containing  multiple  channels,  producing  a  combined  data  rate  in  the  terabits  per  second 
(1012  b/s)  range. 

A  typical  N-channel  WDM/WDMA  system  is  shown  in  Figure  1 .  In  this  figure  a 
WDM  star  network  is  depicted.  It  consists  of  fixed  or  tunable  single-longitudinal-mode 
semiconductor  laser  diodes  (such  as  the  distributed  feedback  (DFB)  laser  for  the 


Figure  1 .  Dense  WDM  network  model  [From  Ref.  1]. 


transmitters),  a  star  coupler,  optical  amplifiers,  and  optical  filters  (such  as  fixed  tuned  or 
tunable  Fabry-Perot  (FP)  filters)  at  the  receiver  ends  [Ref.  1],  In  this  configuration,  each 
transmitter  produces  a  signal  on  a  different  wavelength  (X,,  X^, X^).  The  signals  are 
then  passed  through  the  star  coupler,  where  the  signals  are  combined  and  distributed  to 
the  receivers.  At  each  of  the  receiver  branches,  an  optical  amplifier  amplifies  the  combined 
wavelengths  (X„  X2, X^).  The  combined  signal  is  then  passed  to  a  direct-detection 
receiver.  The  direct-detection  receiver  uses  an  optical  filter  to  pass  the  desired 
wavelength,  isolating  the  wavelength  of  interest  (Xa,  X^  ....,  X^.  In  this  configuration  any 
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of  the  source  transmitters  can  talk  to  any  receiver  by  tuning  to  the  receiver's  preasssigned 
wavelength. 

Both  ON-OFF  keying  (OOK)  and  frequency  shift-keying  (FSK)  are  feasible  for 
WDMAVDMA  systems.  Direct  current  OOK  modulation  of  the  laser  transmitter 
introduces  frequency  excursions  in  the  laser  frequency  or  "chirp".  When  passed  through 
the  optical  filter,  the  chirped  signal  experiences  significant  degradation  when  the  filter 
bandwidth  is  smaller  than  the  chirp  signal.  Figure  2  depicts  the  filter  responses  of  a 
multi-channel  system.  In  the  ideal  system,  the  optical  filters  isolate  the  signal  spectra  in 
each  channel  (Figure  2(a)).  However,  when  the  signal  is  broadened  by  the  chirp, 
interference  from  adjacent  channels  can  occur,  as  in  Figure  2(b).  This  interference  is  called 
crosstalk.  It  is  therefore  important  to  be  able  to  predict  the  amount  of  frequency 
broadening  produced  by  the  chirp  effect.  Later  efforts  will  relate  the  frequency  broadening 
to  the  crosstalk  errors. 
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Using  MATLAB  and  SIMULINK,  this  thesis  implemented  a  model  of  a  DFB  laser 
with  chirp.  The  characteristics  of  chirp  were  explored  by  applying  different  drive  currents 
to  the  model.  The  laser  chirp  model  is  described  in  the  next  chapter,  including  diagrams 
explaining  its  implementation  in  MATLAB  and  SIMULINK.  The  procedures  used  in 
collecting  data  are  discussed  in  Chapter  m,  followed  by  an  analysis  of  the  data  in  Chapter 
IV.  Finally,  Chapter  V  reviews  the  MATLAB/SIMULINK  model  and  summarizes  the 
findings  found  in  Chapter  IV. 
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EL  LASER  CHIRP  SIMULATION 


The  chirp  model  is  broken  into  two  parts,  modeling  the  laser  drive  current  and 
modeling  the  interactions  between  the  charged  carriers  and  photons  in  a  single-mode 
distributed  feedback  laser.  Both  the  mathematical  expressions  for  the  model  and  its 
implementation  in  MATLAB  are  discussed  in  this  chapter. 

A.  MODEL  OF  LASER  CHIRP 

The  expressions  used  to  model  the  laser  drive  current,  as  well  as  the  laser  rate 
equations  were  provided  from  Ref.  1 . 

1.  Laser  Drive  Current 

The  current  pulse  shape,  lp(t),  generated  from  the  laser  driver  is  dependent  on  the 
present  bit  and  the  previous  bit.  The  drive  current  is  characterized  by  Equation  (1)  below: 


hit)  = 


Ibias  +  /m[l-exp(-^) 

I  bias  +  /m[exp  (  Xr” ^ 


1  bias 

1  bias  +1  m 


if  present  bit  =  1,  previous  bit  =  0 

if  present  bit  =  0,  previous  bit  =  1 

if  present  bit  =  0,  previous  bit  =  0 
if  present  bit  =  1,  previous  bit  =  1 


(1) 


where  lbias  is  the  dc  bias  current,  Im  is  the  modulation  current,  and  rr  is  the  rise  time.  The 
fall  time  is  assumed  to  be  equal  to  the  rise  time.  Figure  3  shows  samples  of  different 
waveforms  created  using  Equation  (1).  The  waveform  at  the  top  of  Figure  3  has  an 
alternating  bit-pattern  of  ONEs  and  ZEROs.  When  the  first  ONE  is  the  "present  bit",  the 
waveform  rises  from  Ibias  to  7m.  When  the  next  bit,  a  ZERO  bit,  becomes  the  "present  bit", 
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the  waveform  starts  to  descend  toward  Ibias.  The  current  continues  to  ascend  and  descend 
in  this  manner  for  all  of  the  remaining  bits  in  the  bit-pattern.  In  the  bottom  waveform  of 
Figure  3,  the  bit-pattern  consists  of  alternating  pairs  of  ONEs  and  ZEROs.  In  the  first  pair, 
the  first  ONE  bit  causes  the  waveform  to  ascend  as  before.  However,  since  the  second  bit 
is  also  a  ONE,  the  drive  current  remains  at  Im  for  the  second  bit  period.  Similarly,  the  first 
ZERO  bit  causes  the  current  to  descend  toward  Ibms,  while  the  second  ZERO  bit  causes 
the  drive  current  to  hold  at  Ibuts.  The  final  two  ONE  bits  cause  the  current  to  ascend  and 
hold,  repeating  the  actions  caused  by  the  first  pair  of  ONE  bits.  Also  note  the  different 
horizontal  scales  in  Figure  3.  The  shape  of  the  drive  current  is  dependent  on  the  bit  rate 
and  the  scales  on  the  horizontal  axes  differ  to  reflect  the  bit  rates  used. 
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2.  Single-mode  Distributed  Feedback  (DFB)  Laser 


The  laser  rate  equations  [Refs.  2  and  3]  were  used  to  model  the  operation  of  the 
DFB  laser.  These  coupled,  nonlinear  differential  equations  describe  the  interactions 
between  the  photon  density  p(t),  the  carrier  density  n(t),  and  the  phase  of  the  electric  field 

¥t)- 


^-  =  YG{i)[rit)-n0)pit) 
at 


pjj) ,  yrrtfi 

Xp  +  x„ 


(2) 


^r  =  ^fr- °(0  WO  -  no]  Pit)  -  ^ 


(3) 


dm 

dt 


(4) 


Here,  T  is  the  mode  confinement  factor,  no  is  the  earner  density  at  transparency,  rp  is  the 
photon  lifetime,  p  is  the  fraction  of  spontaneous  emission  coupled  into  the  lasing  mode, 
t  is  the  carrier  recombination  lifetime,  q  is  the  electron  charge ,  V a  is  the  active  layer 

n 

volume,  vs  is  the  group  velocity,  ac  is  the  gain  coefficient,  and  a  is  the  linewidth 
enhancement  factor. 

Table  1  lists  the  parameter  values  used  in  the  laser  rate  equations.  The  values 
shown  are  the  same  values  presented  in  Ref.  3.  Table  1  also  shows  the  labels  used  in  the 

computer  simulations  for  future  reference. 

The  parameter  G(t)  is  the  saturable  gain  coefficient,  defined  by 


G(t)  = 


Vg  f£ 
\+ep(t) 


(5) 
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1  Notation 

Definition 

Value 

SIMULINK  label  || 

r 

mode  confinement  factor 

0.4 

Gamma 

electron  density  at  transparency 

1018  cm'1 

n_0 

rp 

photon  lifetime 

3  ps 

tau_p 

*■„ 

electron  lifetime 

1  ns 

tau_n 

p 

spontaneous  emission  factor 

3  x  10 5 

beta 

9 

electron  charge 

1.6  x  10-'9  C 

q 

K 

active  volume 

1.5  x  10'10cm3 

V_a 

a 

linewidth  enhancement  factor 

5 

alpha 

v* 

group  velocity 

8.5  xlO9  cm/s 

v_9 

oc 

gain  coefficient 

2.5  x  10'“  cm2 

a_0 

total  quantum  efficiency 

0.4 

eta_0 

e 

gain  suppression  factor 

1017  cm3 

eps 

Table  1.  Definition  and  values  of  model/simulation  parameters  [Ref.  3], 


where  e  is  the  gain  compression  factor.  The  value  of  e  is  also  given  in  Table  1 . 

To  find  the  optical  power  emitted  and  instantaneous  chirp  frequency,  Equations 
(2)  -  (5)  are  numerically  integrated  using  the  drive  current  produced  from  Equation  (1). 
The  optical  power  P(t)  is  related  to  the  photon  density  p(t)  by 


n/j\  Pi1)  Var\ohv 
n  )~  2  Tzp 


(6) 


Similarly,  the  instantaneous  chirp  frequency,  Av(t),  is  found  from 


Av(r)  = 


i  dm 

2n  dt 


(7) 


where  r\0  is  the  differential  quantum  efficiency  of  the  laser.  The  value  for  rjo  is  also  given  in 
Table  1. 
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Finally,  substituting  Equations  (2)  and  (4)  into  Equation  (7),  the  frequency  chirp, 
Av(t),  can  be  expressed  as: 


a 

[  i  <m 

4jt 

Pif)  dt 

piXOT 

tn/KO  J. 


(8) 


The  "transient  chirp",  represented  by  the  first  term  in  Equation  (8),  occurs  when 
the  laser  is  turned  off  and  on.  The  sudden  change  in  current  density  disrupts  the 
equilibrium  of  the  system,  creating  a  series  of  relaxation  oscillations  on  the  carrier  density 
signal.  This  ultimately  affects  the  refractive  index  of  the  lasing  medium  and  causes  the 
laser  emission  wavelength  to  shift  correspondingly. 

The  "adiabatic  chirp",  represented  by  the  second  term  in  Equation  (8),  is  the  steady 
state  frequency  offset  between  on  and  off  levels  which  exist  because  of  differences  in 
carrier  density  and  gain  compression  levels.  [Ref.  1] 

Several  studies  have  shown  the  ability  to  reduce  the  frequency  chirp.  For  example, 
the  linewidth  enhancement  factor,  a,  can  be  reduced  from  approximately  a  value  of  5  (in  a 
typical  DFB  laser)  to  approximately  3,  either  by  detuning  the  emission  wavelength  [Refs. 
4-6]  or  by  using  a  multiple-quantum  well  (MQW)  structure  [Refs.  7  and  8],  Additionally, 
other  parameters  such  as  T  and  Va  can  have  a  lesser  effect  on  the  size  of  the  frequency 
chirp  [Ref.  1], 


B.  MATLAB  SIMULATION  OF  CHIRP  MODEL 

MATLAB  and  SIMULINK  were  used  to  simulate  the  model.  A  description  of 
SIMULINK  is  given  later.  MATLAB  is  a  computer  program  that  performs  mathematical 
computations.  During  a  MATLAB  session,  the  user  performs  calculations  by  executing 
instructions  in  a  MATLAB  workspace.  The  results  of  the  calculations  are  stored  within 
the  workspace.  One  of  the  features  of  MATLAB  is  the  ability  to  save  a  workspace  in 
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computer  memory.  This  allows  the  user  to  recall  the  workspace  and  restart  the  MATLAB 
session  using  the  same  variables. 

The  user  can  also  automate  MATLAB  by  using  "m-files".  An  m-file  is  a  set  of 
ordered  MATLAB  instructions,  including  calls  to  other  m-files.  An  m-file  can  either  be 
executed  to  perform  a  single  task  or  be  defined  as  a  function.  When  defined  as  a  function, 
an  m-file  can  be  used  like  any  other  MATLAB  function,  operating  on  the  inputs  provided 
to  it. 

The  following  sections  describe  the  computer  implementation  of  the  laser  chirp 
model.  The  laser  drive  current  was  created  using  MATLAB  m-files  and  is  described  first. 
Then,  after  a  description  of  SIMULINK,  the  simulation  model  of  the  DFB  laser  is  outlined 
and  discussed. 

1.  MATLAB  Implementation  of  the  Laser  Drive  Current 

In  Equation  (1),  the  laser  drive  current  for  a  steady  sequence  of  ZERO  bits  was 
I  bias-  Likewise  the  drive  current  for  a  steady  sequence  of  ONE  bits  was  Ibias  +  Im.  In  the 
MATLAB  implementation  of  Equation  (1),  the  drive  current  for  a  steady  sequence  of 
ONE  bits  was  labeled  as  I  on.  The  value  of  I  on  was  defined  as  1.85  times  the  threshold 
current.  (The  threshold  current,  1  thresh,  is  the  drive  current  required  for  the  diode  to 
begin  lasing.  In  this  thesis,  I  thresh  is  assumed  to  have  the  value  of  33.5  mA.)  Similarly,  in 
the  MATLAB  implementation  of  Equation  (1),  the  drive  current  for  a  steady  sequence  of 
ZERO  bits  was  labeled  as  I_off 

l_off  varied  in  the  computer  simulations.  Its  value  was  a  function  of  the  threshold 
current  and  the  "extinction  ratio".  The  extinction  ratio,  r,  is  defined  as  the  ratio  of  lower 
optical  power  to  higher  optical  power  in  a  digitally  modulated  scheme  [Ref.  9],  The 
following  discussion  shows  the  relationship  of  the  MATLAB  implementation  of  Equation 
(1),  with  the  expressions  given  earlier.  As  shown  below,  the  MATLAB  implementation  of 
Equation  (1)  involved  the  incorporation  of  the  extinction  ratio  and  threshold  current. 
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Equation  (9)  shows  the  extinction  ratio  defined  as  a  function  of  I_off,  I_on,  and 
7  thresh. 


Power  off 
PoweroN 


Ijm  -  Ijhresh 
-  l_off  -  Ijhresh 
0 


for  I_off£  I_thresh 
for  I_ofF  <  I_thresh 


(9) 


Solving  Equation  (9)  for  I  off  and  substituting  I  on  with  1 .85*/  thresh  results  in 
Equation  (10), 


I  jiff  =  Ithresh{\  +  0.85r).  (10) 

The  following  equations  relates  Ibjas  and  Im  with  I  jiff,  I_on,  I  thresh  and  r.  First,  we 
know  that 


I  bias  -  I_off  =  Ijhresh{\  +  0.85r). 


OD 


Since 


I  on  =  I  bias  +Im  =  1  -85  *  I  thresh , 


(12) 


Im  can  be  put  in  terms  of  I  thresh  and  r  by  solving  for  Im  and  substituting  for  I  on  and  Ibias, 

Im  =  I  thresh  *  .85(1  -  r).  (13) 


Finally,  substituting  for  the  numerical  value  of  I  thresh  results  in  Equation  (14)  and  (1 5). 
These  expressions  relate  the  Ibias.  and  Im  to  the  variables  used  in  the  MATLAB 
implementation  of  Equation  (1), 


Ibias  ^  33.5(1  +.85r)mA 


(14) 
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and 


Im  =  28.5(1  -r)  mA. 


(15) 


The  MATLAB  implementation  of  the  laser  drive  current  was  produced  in  the  main 
simulation  m-file,  chrp_mn.m  (see  Appendix  A).  Created  as  a  data  array,  the  laser  drive 
current,  1 _p(t),  was  initialized  to  hold  1000  elements.  Based  on  the  length  of  the 
bit-pattern,  the  elements  were  divided  into  three  equal  increments.  As  shown  in  Figure  4, 
these  increments  corresponded  to  the  bit  period,  Tb,  of  the  bit-pattern.  The  total 
simulation,  therefore,  was  the  number  of  bits  in  the  bit-pattern  times  the  bit  rate.  For  each 
bit  in  the  bit-pattern,  chrp  mn.m  determined  the  value  of  the  present  bit  and  the  previous 
bit  and  applied  the  MATLAB  implementation  of  Equation  (1).  The  values  for  the 
waveform  were  calculated  for  the  entire  bit  period  (increment)  and,  as  shown  in  Figure  4, 
each  increment  was  appended  to  the  previous  increment.  To  start  the  creation  of  the  1 _p(t) 
data  array,  an  initial  bit,  init,  was  required.  Init  simulates  the  "previous  bit"  for  the  first  bit 
in  the  bit-pattern. 


Figure  4.  Generation  of  the  laser  drive  current. 
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Examples  of  different  drive  currents  are  plotted  in  Figure  5.  Figure  5  shows  three 
waveforms,  with  different  bit  rates,  bit-patterns  and  extinction  ratios.  As  in  Figure  3,  the 
shape  of  the  drive  currents  ascend  or  descend  based  on  the  bit-pattern.  Additionally,  the 
scales  on  the  horizontal  axes  reflect  the  use  of  different  data  rates.  Figure  5,  however, 
shows  the  effect  of  the  extinction  ratio  on  the  shape  of  the  laser  drive  current.  By  holding 
I  thresh  (and  therefore  I  on)  constant,  the  value  of  I  off  increases  when  the  extinction 
ratio  increases.  Although  the  bottom  plot  shows  an  extinction  ratio  of  0.8235  for 
illustration  purposes,  typical  values  for  OOK  modulation  would  be  lower.  For  OOK 
modulation,  Ibms  is  set  close  to  I  thresh  (i.e.,  at  a  low  extinction  ratio),  resulting  in  higher 
signal-to-noise  ratios  or,  equivalently,  lower  bit  error  rates. 


Bit-pattern  =>  1  001 


2  4  6 

time  (ns) 


Mhresh  =>  33.5  mA 
|_on  =>  61 .98  mA 

I _ off  =>  36.85  mA 

Extinction  Ratio  =>  0.1 176 
Bit  rate  =>  0.5  GHz 


Bit-pattern  =>  1  010 


0.4  0.6 

time  (ns) 


Mhresh  =>  33.5  mA 
l_on  =>  61 .98  mA 
l_oH  =>  43.55  mA 
Extinction  Ratio  =>  0.3529 
Bit  rate  =>  3  GHz 


Bit-pattern  =>0110 


0.2  0.3 

time  (ns) 


Mhresh  =>  33.5  mA 
l_on  =>  61 .98  mA 
l_otf  =>  56.95  mA 
Extinction  Ratio  =>  0.8235 
Bit  rate  =>  6  GHz 


Figure  5.  Examples  of  MATLAB  generated  laser  drive  currents. 


After  I jp(t)  was  created,  chrp_mn.m  passed  the  drive  current  data  array,  as  well  as 
the  corresponding  time  data  array,  to  the  SIMULINK  DFB  laser  model,  chirp,  m.  The  next 
section  provides  a  brief  overview  of  SIMULINK  and  outlines  the  creation  of  the  laser 
model. 


2.  Single-mode  DFB  Laser  SIMULINK  Simulation 

SIMULINK  is  a  MATLAB  tool,  designed  to  transform  a  graphical  representation 
of  a  mathematical  model  into  a  MATLAB  m-file.  The  SIMULINK  process  has  two  steps. 
First,  the  model  is  graphically  defined  and,  then,  the  model  is  simulated  (analyzed). 

Like  other  SIMULINK  simulations,  the  DFB  laser  model  was  defined  graphically 
by  connecting  "function  blocks"  together.  Depicted  in  Figure  6,  a  "function  block" 
performs  a  mathematical  operation  on  its  inputs,  a  and  b ,  to  produce  an  output,  ffa,  b). 
Several  function  blocks  were  provided  in  SIMULINK.  The  ones  most  commonly  used  for 
this  model  were  sum  blocks,  product  blocks  and  inverting  blocks.  Sum  blocks  were  used 
to  perform  addition  or  subtraction  on  their  inputs.  Product  blocks  multiplied  their  inputs 
together,  and  inverting  blocks  inverted  their  inputs.  The  inverting  blocks  and 
multiplication  blocks  were  used  to  perform  the  division  operation. 


r - ^ 


EJH 

input 

H — * 

f(u.v) 

- +  f(a.b) 

output 

input 

function  block 

example 

inputs  function 

output 

a  =  a 
b  =  2 

addition 

5 

J 

Figure  6.  Example  of  a  SIMULINK  function  block. 
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Other  pre-defined  blocks  were  also  used  in  this  model.  Gain  blocks,  graphically 
represented  by  an  op  amp  symbol,  were  used  to  multiply  an  input  by  a  fixed  amount. 
Blocks  called  "to  workspace"  and  "from  workspace"  were  used  to  pass  data  to  and  from 
a  MATLAB  workspace.  Input  blocks  were  used  to  provide  constant  values  to  the  various 
function  blocks.  Finally,  an  integrator  block,  graphically  depicted  as  a  block  with  a  "1/s" 
label,  performed  numerical  integration. 

Once  the  model  was  graphically  defined,  it  could  be  simulated.  However,  in  order 
to  perform  the  required  numerical  integration,  an  integration  type  and  initial  values  for  the 
integration  needed  to  be  selected.  Several  integration  types  were  provided  in  SIMULINK: 
Runge-Kutta  third  and  fifth  order  methods,  Gear's  predictor-corrector  method,  Adam's 
predictor-corrector  method,  and  Euler's  method.  The  Runge-Kutta  third  order  method 
was  used. 

The  initial  conditions  used  to  integrate  the  derivative  of  the  photon  density  and 
carrier  density  were  also  required.  Two  sets  of  values  were  necessary,  corresponding  to 
two  initial  bit  conditions  (init  =  0  and  init  =  1).  To  find  the  initial  values  when  init  =  0,  an 
all  ZERO  bit  pattern  was  used.  Since  there  is  no  change  in  the  bit-pattern,  the  photon 
density  and  carrier  density  curves  should  remain  at  equilibrium.  However,  if  the  initial 
conditions  for  the  photon  density  and  carrier  density  were  not  set  correctly,  the  signals 
would  fluctuate  until  they  reached  equilibrium.  In  order  to  run  the  simulation,  arbitrary 
initial  conditions  were  provided  for  the  photon  and  carrier  density.  After  the  simulation 
was  completed,  plots  of  the  carrier  and  photon  density  signals  were  examined.  Once  it  was 
determined  that  equilibrium  levels  for  both  photon  density  and  carrier  density  had  been 
reached,  the  equilibrium  values  were  noted  and  used  as  the  initial  condition  values  ( p(0 ) 
and  n(0)).  Likewise  to  find  the  initial  conditions  when  init  =  1,  the  model  was  simulated 
using  an  all  ONE  bit  bit-pattern.  When  equilibrium  levels  had  been  reached,  the  last  value 
of  the  photon  density  and  carrier  density  were  noted  and  used  as  the  initial  values,  p(l) 
and  n(l).  Table  2  shows  the  initial  conditions  that  were  determined  and  used  in  the 
simulations. 
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init  =  0 

init  =  1 

p(0)  =  .271256el5 
n(0)  =  1.39315el8 

p(l)=  1.42154el5 
n(l)=  1.39772el8 

photon  density 
carrier  density 


Table  2.  Initial  conditions  used  in  simulation. 


One  of  the  features  of  SIMULINK  was  the  ability  to  group  blocks  together  to 
form  "subsystems".  Subsystems  divide  the  model  into  small  logical  groups.  Once  a 
subsystem  was  defined,  it  was  connected  in  the  same  manner  as  the  pre-defined  function 
blocks.  Inside  a  subsystem  block,  were  "inport"  and  "outport"  blocks.  These  blocks 
connect  inputs  and  outputs  of  the  subsystem  to  the  blocks  within  the  subsystem. 

Figure  7  shows  the  SIMULINK  graphical  model  of  the  laser,  which  consisted  of 
five  major  subsystems:  the  saturable  gain  coefficient  (G  (t ) ),  the  derivative  of  the  photon 
density  (dn/dt),  the  derivative  of  the  carrier  density  (dp/dt),  the  instantaneous  chirp 
frequency  (Chirp),  and  the  optical  power  emitted  (Power). 


Figure  7.  Major  sections  of  simulation  model. 


In  reviewing  the  block  diagrams,  the  arrows  indicate  the  direction  of  data  flow. 

For  clarity,  all  the  inputs  to  a  subsystem,  or  a  function  block,  are  shown  on  the  left  and  the 
outputs  on  the  right.  Constants  used  in  the  simulation,  found  in  Table  1,  are  indicated  by 
their  respective  values  inside  an  input  block  and  are  further  distinguished  by  a  shadow 
effect. 

In  order  to  help  follow  the  SIMULINK  diagrams,  a  detailed  description  of  the 
saturable  gain  coefficient  subsystem  (Figure  8)  is  given.  The  other  major  subsystems  are 
created  similarly  and  are  shown  in  Figures  9-12. 


a.  Saturable  Gain  Coefficient,  G(t) 

Figure  8  shows  the  breakdown  of  the  saturable  gain  coefficient  defined  by 
Equation  (5).  This  equation  is  repeated  in  the  upper  left  comer  of  the  figure.  The  diagram 
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at  the  bottom  left  of  the  figure  represents  the  overall  subsystem  G  (t) .  This  diagram 
shows  how  the  overall  subsystem  is  further  separated  into  two  "sub  subsystems"  and  a 
product  block.  From  Equation  (5),  G(t)  was  found  by  taking  the  product  of  a0  and  vg,  and 
multiplying  by  the  inverse  of  1  +  sp(t).  The  two  multiplicands  are  produced  in  the  two 
"sub  subsystems",  shown  at  the  right  side  of  Figure  8. 

The  subsystems  for  G  ( t )  are  labeled  a_0  *  v_g  and  1  /  ( 1 + ep  s  *  p  ( t). 
The  formation  of  the  two  subsystems  are  now  considered. 

1.  The  top  right  subsystem,  a_0*  v_g,  calculates  the  product  of  two 
constants  a_0  and  v_g,  and  the  result  is  fed  into  the  product  block  of 
the  overall  diagram. 

2.  The  bottom  right  figure  is  similarly  broken  up: 

a)  The  constant  eps  and  the  input  variable  p  ( t )  are  multiplied 
together  by  a  product  block. 

b)  The  result,  eps*p  ( t),  is  then  added  to  a  constant,  1,  in  a 
sum  block. 

c)  The  result  of  the  addition  is  then  passed  through  an  inverting 
block,  yielding  1/ (l+eps*p  (t)  ) . 

d)  Finally,  l/(l+eps*p(t))  is  fed  into  the  product  block  of 
the  overall  diagram. 

b.  Derivative  of  Carrier  Density,  dn(t)/dt 

Similarly,  Figure  9  shows  the  connections  that  created  the  derivative  of  the 
carrier  density.  The  overall  subsystem,  displayed  at  the  left  side  of  Figure  (9),  shows  three 
subsystems  and  a  product  block  feeding  into  a  sum  block.  From  Figure  7,  four  input 
variables  are  supplied  to  the  dn/dt  overall  subsystem:  G  ( t),  p  ( t),  n  ( t )  and 
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I  p  ( t).  These  input  variables  are  further  distributed  to  the  subsystems  within  the  dn/  dt 
subsystem. 


Figure  9.  Derivative  of  the  carrier  density  subsystem. 

Equation  (3),  repeated  at  the  top  of  Figure  9,  was  used  to  create  the 
dn/dt  subsystem.  The  first  term  on  the  right  side  of  Equation  (3)  was  calculated  in  the 
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top  right  subsystem  diagram.  This  subsystem  takes  the  input  variable  I_p  ( t )  from  an 
inport  block  and  multiplies  its  value  with  the  inverted  product  of  q  and  V_a.  The  middle 
term  on  the  right  side  Equation  (3)  was  split  into  two  parts.  The  first  part  was  calculated 
by  the  subsystem  depicted  in  the  middle  right  of  Figure  9.  Here,  the  input  variables  G  ( t ) 
and  p  ( t )  are  multiplied  with  the  constant  n_0  in  a  product  block.  The  second  part  of  the 
middle  term  was  calculated  with  the  product  block  found  in  the  overall  dn/dt  subsystem 
diagram.  The  product  block  multiplied  the  input  variables  G  (t) ,  n  (t )  and  p  ( t ) ,  and 
fed  the  result  directly  into  the  sum  block  of  the  dn/dt  subsystem.  The  last  subsystem, 
shown  at  the  bottom  right  of  Figure  (9),  calculated  the  last  term  of  Equation  (3).  It  took 
the  product  of  the  input  variable  n  ( t )  and  multiplied  it  by  the  inverted  value  of  t  au_n. 
The  output  of  this  subsystem,  as  with  the  other  subsystems  shown  on  the  right  side  of 
Figure  (9),  were  fed  into  the  sum  block  in  the  overall  dn/  dt  subsystem  diagram.  (Note 
that  addition  and  subtraction  was  indicated  by  a "+"  or sign  on  the  sum  block.)  The 
mathematical  operation  performed  in  the  sum  block  reflected  the  same  operation  required 

in  Equation  (3). 

c.  Derivative  of  Photon  Density,  dp(t)/dt 

Figure  10  shows  the  composition  of  the  derivative  of  the  photon  density. 
Similar  to  the  dn/dt  subsystem,  this  subsystem,  however,  used  only  three  input 
variables:  G  ( t ) ,  p  ( t ) ,  and  n  ( t ) .  Using  Equation  (2)  (repeated  at  the  top  of  the  Figure 
10),  the  dp/dt  subsystem  was  further  broken  into  three  subsystems  and  a  product  block. 
The  diagram  on  the  left  side  of  Figure  10  shows  the  overall  dp/dt  subsystem. 

The  first  term  on  the  right  side  of  Equation  (2)  was  split  into  two  parts. 
The  first  part  was  calculated  by  a  product  block,  located  in  the  overall  dp/ dt  subsystem. 
This  product  block  multiplied  Gamma  with  the  input  variables  G  (t) ,  n  (t )  and  p  ( t). 
The  result  of  the  product  block  fed  directly  into  the  sum  block  of  the  overall  dp/  dt 
subsystem.  The  second  part  of  the  first  term  was  calculated  in  the  top  right  subsystem. 
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Figure  10.  Derivative  of  the  photon  density  subsystem. 


Here,  the  input  variables  G  ( t )  and  n  ( t )  are  multiplied  with  the  constants  Gamma  and 
n  0.  The  secondterm  on  the  right  side  of  Equation  (2)  was  calculated  in  the  middle  right 
subsystem.  The  input  variable  p  ( t )  was  multiplied  by  the  inverted  value  of  t  au_p. 
Finally,  the  third  term  of  Equation  (2)  was  calculated  by  the  bottom  right  subsystem.  Here, 
the  input  variable,  n  (t),  was  multiplied  by  beta,  Gamma  and  the  inverted  value  of 
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tau_n.  All  three  subsystems,  shown  at  the  right  of  Figure  10,  were  fed  into  the  sum 
block  of  the  overall  dp/dt  subsystem.  As  per  Equation  (2),  the  sum  block  indicates 
whether  addition  or  subtraction  was  performed. 

After  the  derivative  of  the  photon  density  and  carrier  density  were 
calculated,  their  outputs  were  fed  into  integrators.  The  integrators  calculated  n  ( t )  and 
p  ( t).  As  shown  in  Figure  7,  n  ( t )  and  p  ( t )  are  distributed  and  fed  back  into  the 
system.  The  initial  values  used  for  the  integration  are  given  in  Table  2,  as  discussed  earlier. 

d.  Instantaneous  Chirp  Frequency,  Av(t) 

The  culmination  of  the  SIMULINK  model  was  the  calculation  of  the 
instantaneous  chirp  frequency,  shown  in  Figure  1 1 .  Here  Equations  (7)  and  (10)  were 
combined  to  form 


Av(0  =  £(§)  (rvsa„K0  -».]-£).  (15) 

Equation  (15)  was  used  to  create  the  instantaneous  chirp  frequency  subsystem. 

Starting  at  the  left  side  of  Figure  1 1,  the  constant  n_0  was  subtracted  from 
the  input  variable  n  ( t )  in  a  sum  block.  The  result  of  this  first  sum  block  was  then 
multiplied  with  the  constants  v_g,  a_0  and  Gamma  in  a  product  block.  The  result  of  this 
first  product  block  was  fed  into  a  second  sum  block.  The  second  sum  block  took  the 
difference  between  the  result  of  the  product  block  and  the  inverse  value  of  the  constant 
tau  p.  The  result  of  the  second  sum  block  was  then  multiplied  in  a  second  product  block 
by  a  factor  of  ot/2  (a/2  was  calculated  in  a  gain  block).  The  result  of  the  second  product 
block  was  scaled  by  a  second  gain  block.  The  scaling  factor  for  this  gain  block  was  set  at 
1/(2ti). 
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Figure  1 1 .  Instantaneous  chirp  frequency  subsystem 


e.  Optical  Power  Emitted  P(t) 


The  optical  power  emitted  was  calculated  as  in  Figure  12.  By  substituting 
the  frequency,  v,  for  the  speed  of  light  divided  by  the  wavelength,  v  =  c/X,  Equation  (6) 
becomes 


P,A  PV)  VgT\ohc 

K)  2  TxpX 


(16) 


Starting  at  the  left  side  of  Figure  12,  the  product  of  all  the  denominator 
terms  (lambda,  2*  tau_p,  and  Gamma)  of  Equation  (16)  were  calculated  in  a  product 
block.  The  result  of  the  product  block  was  then  inverted  and  multiplied  with  the 
numerator  terms  (input  variable  p  (t ) ,  V_a,  eta_0,  h,  and  c)  in  a  second  product 
block. 
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Figure  12.  Optical  power  emitted  subsystem. 

The  instantaneous  chirp  and  the  power  emitted  subsystems  were  connected 
to  "to  workspace"  blocks.  As  mentioned  earlier  these  blocks  and  the  "from  workspace" 
block  pass  data  between  the  MATLAB  workspace.  Created  as  data  arrays,  the  emitted 
power  (Power),  the  instantaneous  chirp  frequency  (chrp),  the  carrier  density  (n),  the 
photon  density  (p),  and  a  copy  of  the  laser  drive  current  ( i_p)  (shown  in  Figure  7)  were 
all  saved  in  a  workspace. 

In  this  chapter  the  mathematical  expression  for  the  laser  drive  current,  as 
well  as  the  laser  rate  equations,  were  modeled  into  MATLAB  and  SIMULINK.  The  next 
chapter  describes  the  collection  and  storage  of  data,  and  how  placing  the  data  in  separate 
workspaces  automated  the  process. 
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HL  DATA  COLLECTION 


The  procedure  for  gathering  data  was  to  1)  create  input  parameters  for  the  laser 
drive  current,  2)  create  the  drive  current,  3)  set  the  initial  values  for  integration,  4)  run  a 
simulation  of  the  DFB  laser  using  the  drive  current  as  an  input,  and  5)  collect  and  display 
the  data  produced  by  the  simulations.  As  described  in  Chapter  n,  chrp  mn.m  and  chirp. m 
accomplished  steps  2  and  4.  The  initial  values  in  Step  3  are  given  in  Table  2.  However,  to 
accomplish  the  other  steps,  a  "workspace  structure"  was  used.  In  a  workspace,  the 
simulation  variables  and/or  output  data  arrays  were  maintained  as  individual  data 
packages.  After  being  assigned  a  name,  the  data  packages  were  saved  in  memory  and 
stored  for  future  use.  When  needed,  the  data  packages  were  reloaded  into  MATLAB  and 
the  information  was  extracted.  This  flexibility  made  it  easier  to  retrieve  and  display 
simulation  data. 

A.  DATA  COLLECTION  TO  FIND  CHARACTERISTICS  OF  CHIRP 

The  simulation  was  used  in  two  ways:  first,  to  find  the  general  characteristics  of 
frequency  chirp  and,  second,  to  see  how  the  extinction  ratio  value  affects  the  model.  In 
obtaining  a  generalization  of  chirp  characteristics,  the  m-file  stage l.m  was  created. 

Stage l.m  set  up  simulations  by  combining  different  bit-patterns  and  data  rates.  A  three  bit 
bit-pattern  with  an  initial  condition  (m/7)  of  0  or  1  was  used  for  the  simulations.  Table  3 
shows  the  different  initial  conditions  and  bit-pattern  combinations.  Data  rates  of  0.5,  3,  6, 
and  9  gigabits  per  second  (Gb/s)  were  used.  These  values  reflect  data  rates  associated  with 
present  LAN  data  rates  and  future  designs  using  WDM.  The  sixteen  bit-pattern 
combination  and  the  four  data  rates  resulted  in  sixty-four  cases  being  simulated. 

In  the  first  set  of  simulations  the  extinction  ratio,  r,  was  maintained  as  a  constant. 
An  extinction  ratio  of  0.1882  was  assumed.  Since  the  extinction  ratio  remained  the  same, 
Ioff  was  also  constant.  The  resulting  value  for  I  off  was  38.36  mA.  As  mentioned  earlier, 
I  on  was  held  constant  throughout  all  the  simulations. 
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Initial  Condition/Bit-patterns  Combinations  J 

initial  condition 

drive  current 
bits 

initial  condition 

drive  current 
bits 

0 

000 

1 

000 

0 

001 

1 

001 

0 

010 

1 

010 

0 

Oil 

1 

Oil 

0 

100 

1 

100 

0 

101 

l 

101 

0 

110 

1 

110 

0 

in 

1 

ill 

Table  3.  Bit-patterns  used  in  simulation. 


For  each  of  bit-pattem/data  rate  combinations,  stage l.m  simulated  the 
MATLAB/SIMULINK  DFB  laser  model.  The  model  was  initiated  by  a  function  call  to 
chrpmn.m.  After  chrpmn.m  completed  the  simulation  run,  all  the  parameters  which 
defined  the  laser  drive  current  ( r ,  l_on,  Ijhresh,  the  bit  rate,  the  bit-pattern,  and 

init)  were  saved  in  a  workspace  along  with  the  data  arrays  produced  in  the  simulation 
(/' _p,  n,  p Power,  chrp ).  The  workspace  was  then  given  a  name  (e.g.,  easel. mat, 
case2.mat,  etc.)  and  saved  to  memory.  (See  Appendix  A.) 

The  data  arrays  from  the  simulations  were  plotted  and  presented  using  chrpdsp.m 
(see  Figure  13  for  an  example).  In  the  figure,  the  laser  power  emitted,  the  frequency  chirp, 
the  drive  current,  the  photon  density  and  the  carrier  density  were  all  plotted  with  respect 
to  time.  Additionally,  case-specific  information  was  included.  The  bit-pattern  and  bit  rate 
are  given  above  the  plot  of  the  power  emitted  and  frequency  chirp.  The  extinction  ratio,  r, 
along  with  the  values  for  Ijhresh,  I_on,  and  l_off  were  provided  in  the  text  section 
(located  right  of  the  drive  current  plot).  The  maximum  and  minimum  frequency  chirp 
excursions  were  also  calculated  from  the  data  and  are  provided  in  the  text  section.  These 
values  represent  the  most  positive  and  most  negative  frequency  excursions.  The  initial 
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Figure  13.  Display  of  simulation  data. 
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Frequecy  Chirp  (GHz)  [-] 


condition  for  the  photon  density  and  carrier  density  were  also  provided  in  the  title  of  their 
plots. 

It  was  useful  to  display  the  laser  output  power  and  the  frequency  chirp  on  the  same 
plot.  This  was  accomplished  with  a  MATLAB  function  obtained  from  Mathworks, 
plotyy.m  (see  Appendix  AV 

Another  m-file,  stage4.m,  was  also  created  to  help  display  the  data.  Stage4.m 
allowed  the  user  to  call  up  a  particular  case,  or  a  range  of  cases,  and  display  the  data  from 
the  respective  simulations  (see  Appendix  AV 

Analysis  of  the  plots  will  be  provided  later.  However,  as  mentioned  earlier,  chirp 
can  cause  system  degradation  and  limit  the  spacing  in  a  multi-channel  system.  Therefore, 
the  second  area  of  research  focused  on  the  extinction  ratio  and  how  it  affects  the 
frequency  excursions. 

B.  DATA  COLLECTION  FOR  THE  EXTINCTION  RATIO  EFFECTS 

The  extinction  ratio  was  varied  to  examine  its  affect  on  the  chirp.  This  was  done 
by  creating  the  m-file  extnjtr.m  (see  Appendix  AV  Similar  to  stagel.m,  extn  itr.m  was 
used  to  set  up  the  MATLAB/SIMULINK  simulation.  The  different  simulations  involved 
varying  the  extinction  ratio  for  a  given  bit-pattern  and  data  rate.  The  bit-patterns  were 
limited  to  0  1 1 1,  0  101,  1  010  and  1  000.  The  four  bit-patterns  were  combined  with  the 
four  data  rates  previously  used.  The  extinction  ratio  was  varied  between  0  and  1 .  The  step 
size  used  to  vary  the  extinction  ratio  was  0.01 .  This  resulted  in  a  total  of  1,600 
simulations. 

Instead  of  putting  each  simulation  in  a  separate  case,  the  maximum  frequency 
excursions  from  each  of  the  simulations  were  stored  along  with  their  respective  extinction 
ratio  value.  These  values  were  grouped  together  based  on  the  bit-pattern  and  data  rate. 
Each  data  array  held  values  for  simulations  having  the  same  data  rate  and  bit-pattern.  A 
total  of  16  extinction  ratio  data  arrays  were  produced. 
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The  data  from  the  second  group  of  simulations  was  used  to  plot  the  maximum 
frequency  excursion  versus  extinction  ratio  curves.  Figure  14  shows  a  typical  plot  of  the 
maximum  frequency  chirp  versus  extinction  ratio.  In  this  particular  case,  the  plot  shows 
that  the  maximum  frequency  excursions  curve  had  a  minimum  value.  This  minimum  value 
will  be  discussed  later  in  more  detail. 

To  plot  the  extinction  ratio  curves,  a  second  plotting  routine,  plt_yec.m,  was 
created.  Plt  vec.m  used  basic  MATLAB  plotting  features  to  manipulate  and  display  the 
data  from  the  extinction  ratio  data  arrays  (see  Appendix  A). 


Figure  14.  Typical  frequency  versus  extinction  ratio  curve. 


In  this  chapter  two  MATLAB  m-files  were  described,  stage  I  .m  and  extnjtr.m. 
These  m-files  set  up  simulations  of  the  MATLAB/SIMULINK  model.  The  data  that  was 
collected  from  the  model  were  best  analyzed  in  their  graphical  formats.  The  next  chapter 
analyzes  these  plots  in  order  to  make  general  observations  of  the  chirp  effect. 
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IV.  ANALYSIS  OF  DATA 


The  next  two  sections  discuss  the  conclusions  reached  by  analyzing  the  plots 
produced  by  the  simulations.  The  first  set  of  simulations  was  analyzed  in  an  attempt  to  find 
general  characteristics  of  chirp.  These  simulations  focused  on  the  trends  in  the  frequency 
excursions  caused  by  different  bit-patterns  and  data  rates.  The  second  set  of  simulations 
was  designed  to  evaluate  the  effects  of  the  extinction  ratio.  Specifically,  the  effects  of  the 
extinction  ratio  on  the  maximum  frequency  excursions  were  explored. 

A.  BIT-PATTERN/DATA  RATE  ANALYSIS 

The  plots  produced  (provided  in  Appendix  B)  and  the  largest  frequency  excursion 
for  each  case  (listed  in  Table  4)  were  examined  to  find  trends  in  the  laser  chirp.  Although 
some  findings  were  not  consistent  throughout  all  of  the  simulations,  some  general 
observations  could  be  made  on  the  characteristics  of  the  chirp. 

1.  General  Characteristics  of  the  Chirp 

The  differences  in  the  maximum  frequency  excursions  caused  by  a  ZERO  to  ONE 
transition  and  a  ONE  to  ZERO  transition  was  first  noted.  The  plots  revealed  that  the 
transition  from  ZERO  to  ONE  caused  a  greater  frequency  excursion  than  the  transition 
from  ONE  to  ZERO.  In  Figure  15,  the  frequency  chirp  curve  is  plotted  for  a  bit-pattern 
of  1  010  with  a  data  rate  of  0.5  Gb/s.  In  this  bit-pattern,  both  the  transitions  from  ZERO 
to  ONE  and  ONE  to  ZERO  occur.  Here,  the  maximum  frequency  excursion  occurs  when 
the  second  bit  is  processed,  or  more  specifically,  during  the  ZERO  to  ONE  transition. 

Also  note  that  the  transient  chirp  related  to  the  ZERO  to  ONE  transition  spans  between 
- 1  GHz  to  5  GHz.  The  transient  chirp  for  the  ONE  to  ZERO  transition  was  much  smaller 

and  spanned  less  than  2  GHz. 
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Data  Rate  Gb/s 


Pattern 

0.5 

3 

6 

9 

0  000 

0.34 

0.34 

0.34 

0.34 

0  001 

4.55 

11.01 

13.42 

14.45 

0  010 

4.54 

11.01 

13.42 

1445 

0011 

4.54 

11.01 

13.42 

14.45 

0  100 

4.56 

11.01 

13.42 

14.45 

0  101 

4.56 

11.26 

21.95 

>19.87 

0  110 

4.56 

11.01 

13.42 

14.45 

0  111 

4.56 

11.01 

13.42 

14.45 

1  000 

1.88 

-4.55 

-7.26 

-8. 44 

1  001 

4.54 

12.73 

18.93 

25.1 

1  010 

4.43 

11.55 

20.5 

18.78 

1  011 

4.43 

11.55 

20.5 

18.96 

1  100 

1.88 

-4.54 

-7.26 

-8.44 

1  101 

4.43 

11.55 

20.49 

18.57 

1  110 

1.88 

-4.54 

-7.26 

-8.44 

1  111 

1.88 

1.88 

1.88 

1.88 

Table  4.  Largest  frequency  excursion  in  GHz  for  each  case. 
Negative  values  indicate  that  the  frequency  excursion  had  the 
largest  value  in  the  negative  direction. 


From  the  laser  rate  equations,  the  transient  chirp  and,  therefore,  the  maximum 
frequency  excursion,  was  a  function  of  the  model  parameters  (listed  in  Table  1).  Later 
studies  can  apply  different  parameters  to  simulate  the  various  laser  sources  available.  By 
doing  so,  the  transient  chirp  and  its  effects  on  the  maximum  frequency  excursion  can  be 
predicted  for  specific  optical  systems. 
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Figure  15.  Example  of  bit  transitions.  1 thresh  =  33.5  mA,  I  on  =  62.0  mA, 
l_off=  3X.9  mA  with  an  extinction  ratio  of  0.1  X8. 


2.  Effects  of  Increased  Data  Rates  on  Chirp 

The  most  significant  finding  from  the  first  set  of  simulations  was  the  increased 
maximum  frequency  excursion  caused  by  an  increase  in  the  data  rate.  In  general  this 
observation  held  true;  however,  for  the  1  010  and  1  01 1  bit-patterns,  the  maximum 
frequency  excursion  for  the  9  Gb/s  cases  were  lower  than  the  6  Gb/s  cases. 

Additionally,  in  the  higher  data  rate  cases,  the  bit  would  change  well  before  the 
chirp  oscillations  had  time  to  settle.  As  an  example,  Figure  16  shows  the  same  bit-pattern, 
but  with  different  data  rates.  Keeping  in  mind  that  the  bit  pattern  is  1  010,  notice  that  in 
the  0.5  Gb/s  data  rate  case,  the  chirp  settles  before  a  new  bit  is  processed  and  the  pow'er 
curve  maintains  the  general  shape  of  the  drive  current.  However,  in  the  higher  data  rate 
cases,  the  transient  chirp  does  not  settle  within  the  bit  period  causing  the  power  curves  to 
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fluctuate  throughout  the  bit  period.  Later  efforts  will  determine  the  signal  degradation 
caused  when  the  duration  of  transient  chirp  exceeds  the  bit  period. 


Figure  16.  Chirp  effects  at  different  data  rates,  ljhresh  =  33.5  mA,  I  on  = 
62.0  mA,  l_off=  38.9  mA  with  an  extinction  ratio  of  0.188. 


3.  Predictability  of  the  Effects  of  Chirp 

The  chirp  was  more  predictable  in  the  lower  data  rate  cases  than  in  the  higher  data 
rate  cases.  For  example,  in  the  0.5  Gb/s  cases,  in  all  instances  where  a  ZERO  to  ONE 
transition  occurred,  the  maximum  frequency  excursion  was  approximately  4.5  GHz 
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(±  1  GHz).  However,  as  the  data  rate  increased,  the  deviation  in  the  value  of  the 
maximum  frequency  excursion  also  increased.  In  the  3  Gb/s  cases,  the  maximum  frequency 
excursion  deviated  more  than  1  GHz  in  some  cases  and,  in  the  9  Gb/s  cases,  the  maximum 
frequency  excursion  deviated  as  much  as  S  GHz. 

Additionally,  as  the  bit  pattern  became  more  complex  (i.e.,  more  bit  transitions), 
the  maximum  frequency  excursion  became  less  predictable.  This  is  best  seen  in  the  cases 
for  the  1  010  and  0  101  bit-patterns.  Notice  that  the  values  for  these  cases  (listed  in  Table 
4)  are  markedly  different  than  the  other  values  within  the  same  data  rate  column. 

B.  EFFECTS  OF  EXTINCTION  RATIO 

It  became  evident  that  the  chirp  settling  time  was  dependent  on  the  maximum 
frequency  excursions.  Therefore,  by  minimizing  the  largest  frequency  excursions,  the 
chirp  relaxation  oscillations  would  settle  more  rapidly.  In  the  second  set  of  simulations,  the 
dependency  of  the  chirp  on  the  extinction  ratio  was  explored. 

As  mentioned  earlier  only  four  bit-patterns  were  considered  for  this  portion  of  the 
thesis.  Plotting  the  maximum  frequency  excursion  versus  extinction  ratio  resulted  in  a 
variety  of  curves.  In  Figure  17,  the  data  for  0  111  bit-pattern  was  plotted.  This  diagram 
shows  that  a  minimum  occurred  on  each  curve  for  each  of  the  data  rates  considered.  This 
minimum  was  significant  since  it  provides  a  value  for  the  extinction  ratio  that  minimizes 
the  maximum  frequency  excursions.  The  minimum  occurred  around  the  extinction  ratio 
value  of  0.2.  Figure  17  also  shows,  the  maximum  frequency  excursions  increase  as  the  bit 
rate  increases.  This  observation  was  consistent  with  the  results  found  in  the  first  set  of 
simulations. 

In  the  remaining  extinction  ratio  bit-pattern  cases,  the  plots  (provided  in  Appendix 
C)  did  not  follow  any  specific  trend.  Other  than  reaffirming  the  characteristics  observed 
from  the  first  set  of  simulations,  the  plots  did  not  show  any  predictable  behaviors. 
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However,  it  was  noted  that  the  length  of  the  bit-patterns  does  not  reflect  a  realistic  digital 
system.  Future  study  will  have  to  examine  the  extinction  ratio  using  longer  pseudorandom 
patterns.  A  pseudorandom  bit-pattern  simulates  a  more  realistic  drive  current  and  may 
lead  to  a  better  estimate  of  optimal  extinction  ratio  values. 


Figure  17.  Maximum  frequency  chirp  (Bit-pattern  =>  0111). 


In  analyzing  the  data  generated  by  the  simulations,  it  became  apparent  that  an 
infinite  number  of  cases  could  be  considered.  However,  the  trends  and  observations 
provided  from  the  plots  gave  a  general  representation  of  laser  chirp.  The  next  chapter 
reviews  the  creation  of  the  MATLAB/SIMUL1NK  model  and  summarizes  the  findings  in 
this  chapter. 
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V.  SUMMARY 


This  thesis  implemented  a  model  of  a  distributed  feedback  laser  in  order  to  predict 
the  effects  of  chirp.  The  model,  which  was  based  on  expression  for  a  laser  drive  current 
and  the  laser  rate  equations,  was  implemented  using  MATLAB  and  SIMULINK.  The  laser 
drive  current,  which  represented  a  signal  using  an  OOK  format,  was  modeled  in 
MATLAB.  Coded  in  an  m-file,  the  drive  current  was  shaped  depending  on  the  bit-pattern, 
the  data  rate  and  the  current  parameters  (/  on,  loff,  1  thresh,  init  and  the  extinction 
ratio).  The  drive  current  was  created  as  a  data  array.  Once  the  drive  current  data  array 
was  produced,  it  was  passed  into  the  DFB  laser  SIMULINK  model.  Using  the  drive 
current,  the  model  simulated  the  dynamics  of  the  DFB  laser  producing  data  arrays  for  the 
laser  power  emitted,  the  instantaneous  chirp  frequency,  photon  density  and  carrier  density. 
These  data  arrays,  the  drive  current  data  array  and  parameters  that  shaped  the  drive 
current  were  all  saved  in  a  workspace.  Using  workspaces  allowed  the  data  collection  to  be 
automated.  By  setting  up  the  simulation  in  a  workspace  and  saving  the  data  produced 
from  the  simulation  in  the  same  workspace,  information  was  logically  grouped  in  data 
packages.  After  the  simulations  were  completed,  specific  data  packages  could  be  recalled 
and  the  information  produced  could  be  extracted  and  analyzed. 

Two  simulation  sets  were  conducted.  The  first  simulation  set  was  performed  to 
find  general  characteristics  of  chirp.  These  simulations  focused  on  the  bit  pattern  and  the 
data  rates.  It  was  found  that  the  ZERO  to  ONE  transitions  caused  the  greatest  frequency 
excursions.  This  excursion,  caused  by  the  transient  chirp,  was  a  function  of  the  laser 
parameters  used  (listed  in  Table  1).  Later  studies  can  apply  different  parameters  to  model 
the  various  laser  sources  available.  By  doing  so,  the  transient  chirp  and  its  effects  on  the 
maximum  frequency  excursion  can  be  predicted  better. 

The  most  significant  finding,  however,  was  the  effect  of  increased  data  rates.  As 
the  data  rate  increased,  the  duration  of  the  relaxation  oscillation  in  the  transient  chirp  did 
not  settle  out  before  a  new  bit  was  processed.  Later  efforts  will  need  to  determine  the 
signal  degradation  caused  when  the  duration  of  transient  chirp  exceeds  the  bit  period. 
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Finally,  in  the  first  set  of  simulations,  it  was  determined  that  the  predictability  of 
chirp  decreases  as  the  bit  rate  increases  and  as  the  bit-pattern  becomes  more  complex. 

This  observation  tends  to  emphasis  the  need  to  reduce  the  chirp,  in  order  to  increase  the 
data  rate  of  a  channel. 

The  second  set  of  simulations  provided  some  insight  into  the  effects  of  varying  the 
extinction  ratio.  The  plots  of  the  maximum  frequency  excursions  versus  extinction  ratio 
had  varying  results.  For  the  0  111  bit-pattern  cases,  a  minimum  value  near  r  =  0.2 
occurred  in  each  curve.  This  minimum  represents  a  value  of  the  extinction  ratio  where  the 
maximum  frequency  excursion  would  be  minimized.  The  other  plots,  however,  did  not 
exhibit  any  predictable  behavior.  Future  study  should  incorporate  the  use  of  longer 
pseudorandom  bit-patterns  to  simulate  a  more  realistic  drive  current.  These  simulation 
may  lead  to  a  general  setting  for  the  extinction  ratio  that  would  minimize  the  maximum 
frequency  excursion  for  a  digitally  modulated  signal. 

The  MATLAB  simulation  of  DFB  laser  provided  insight  into  the  behavior  of 
chirp.  The  importance  of  predicting  the  effects  of  frequency  chirp  is  evident  in  the  design 
of  WDMAVDMA  systems  In  a  dense  WDMAVDMA  system,  the  combined  crosstalk  and 
chirp  effect  limit  the  spacing  between  channels  and  degrades  the  transmitted  signal.  By 
characterizing  and  finding  ways  to  limit  the  chirp,  denser  systems  can  be  constructed. 
Combined  with  other  advances  in  optical  technology,  such  as  improved  Fabry-Perot  filters 
and  laser  source  enhancements,  the  high-capacity,  high-data  rate  systems  envisioned  can 
be  realized  and  made  useful  for  mass  data  transfer  applications. 
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APPENDIX  A.  MATLAB  CODE 


This  appendix  contains  MATLAB  codes  used  in  this  thesis.  The  order  of 
presentation  is  given  below: 


A*1 

chrpjnn.m 

main  simulation  program 

E E 

stagel.m 

input  parameter  generator  (different  bit-pattems) 

extn_itr.m 

input  parameter  generator  (different  extinction  ratios) 

A-4 

stage4.m 

display  m-file,  selects  specific  cases 

chrp_dsp.m 

function  call  to  produce  display  graphic 

EB 

vec _plt.m 

display  extinction  ratio  plots 

ES 

plotyy.m 

special  plot  feature,  provided  by  Mathworks 

39 


A-l.  Chrpjmrtrtu 

This  is  the  main  m-file  function  that  initiates  the  simulation.  It  creates  the  laser 
drive  current  and  makes  a  function  call  to  the  SIMULINK  model. 


%chrp_mn.m  This  function  take  input  parameters  from  either 
%  stagel.m  or  extn_itr.m,  and  creates  a  laser  driving 

%  current.  It  then  makes  a  function  call  to  chirp. m  to 

%  execute  a  simulation  run. 

function  [  ]  =  chrp_mn (  ) ; 

global  B  ns  N  a  I_thresh  I_on  I_of f  r  Ijp  Power  chrp  p  . . .  n  i_p  t 
init  case_no; 

%  Calculate  driving  current 
T_B  =  1/B; 
tau_rise  =0.5*T_B; 
tau_fall  =  tau_rise; 

N  =  length (a) ; 

%  Setup  arrays 
ns  =  1000; 

t  =  zeros (ns, 1);  %  Initialize  time  array 

t ( : )  =  linspace ( 0, N*T_B, ns ) ;  %  Divide  N  bit  periods 

%  into  n  pieces. 

I_p  =  zeros (ns, 1);  %  Intialize  I_p  array. 

%%  Start  of  I_p  iteration  routine 

%  Initial  bit 
if  init  ==  1 

I_jstart  =  I_on; 
else 

I_start  =  I_off; 
end 

%  Bit  Pattern 
incr__beg  =  1; 

for  j  =  [IsN] ; 
incr  =  j*ns/N; 

if  a ( 1 , j )  ==  1 
I_aim  =  I_on; 
tau  =  tau_rise; 
else 

I_aim  =  I_off; 
tau  =  tau_fall; 
end 

I_p  (incr_beg:  floor  (incr)  )  =  I__aim  -  (I_aim-I_start )  .*  ... 
exp (-2.2* (t (incr_beg: floor (incr) ) -t (incr_beg) ) /tau) ; 

I_start  =  I__p  (floor  (incr)  )  ;  %save  last  value  of  I_p 
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input  current  parameters 
%%  Bit  period 
%%  Rise-time 
%%  Fall  time 


incr_beg  =  ceil(incr);  %save  position  of  last  I_p 

end 

%%%  End  of  I_p  Creation 


ft%%  Start  SIMULINK  Simulation 

%  Initial  values  for  SIMULINK  Chirp  Model 
start  =  N/B/1000; 
min_inc  =  N/B/1000; 
max_inc  =  N/B/1000; 
tol  =  le-5; 
stop  =  N/B; 

ss  =  [start  stop] ; 

incr  =  [tol  min_inc  max_inc] ; 

%  SYNTAX  NOTE:  xO  must  be  column  array 
if  init  ==  1 

xO  =  [ 1 . 39772* lel8 

1 . 42154*lel5] ; 

else 

xO  =  [ 1 . 39315*lel8 

2 . 71256*lel4] ; 

end 

%  Simulation 

[t,x,y]  =  rk23 (' chirp ' ,  ss,  xO,  incr) ; 

%%%  End  SIMULINK  Simulation 
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A-2.  Stagel.m. 


This  m-file  creates  input  variable  for  chrp_mn.m,  iterating  bit-patterns  and  bit 

rates. 

fcstagel.m  This  m-file  produces  the  input  parameters  for  chrp_mn.m  and 

%  creates  a  workspace  to  store  input  and  output  values. 

clear 

global  B  ns  N  a  I_thresh  I_on  I_off  r  I_p  Power  chrp  p  n  i_j>2  t  ... 

init  case_no; 

fcCalculate  necessary  input  values. 

I_thresh  =  33.5e-3;  %%  I_thresh  is  33.5  mA 

I_on  =  1 . 85*I_thresh;  %%  I_on  is  1.85x  threshold 

x  =  1.16;  I_off  =  x*I_thresh;  %%  Set  I_off 

r= (I_of f-I_thresh) / (I_on-I_thresh) ;  %%  Extinction  ratio. 

A  =  [  0  0  0  %Matrix  with  8  different  patterns 

0  0  1 
0  10 
Oil 
10  0 
10  1 
110 
l  1  1  ]  ; 

b  =  [,5e9  3e9  6e9  9e9] ;  %Four  different  bit  rates 

count  =  0; 
for  i  =  [1:4] 
for  init  =  [0:1] 
for  j  =  [1:8] 

count  =  count  +  1; 
a  =  A  ( j  ,  : )  ; 

B  =  b  (1,  i)  ; 

chrp_mn  %functional  call  to  execute  simulation 

%SETUP  WORKSPACE  AND  STORE  DATA 
case^no  =  count; 

savcase  =  [  ’save  case*  num2str ( count )  1  t  Power  chrp  p  n  i_p2  B  ... 

a  I_thresh  I_on  I_off  r  init  case_no'  ]; 

eval  {[savcase]);  %save  values  in  a  workspace 


end 

end 

end 


%iterate  for  4  Bit  rates 
literate  for  2  initial  conditions 
%iterate  for  8  possible  Bit-patterns 
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A-3.  Extn_itr.  m 

This  m-file  sets  up  the  second  set  of  simulation.  It  iterates  the  extinction  ratio  for 
specific  bit-patterns  and  bit  rates. 

%extn  itr.m  This  m-file  sets  up  the  second  set  of  simulation. 

%  It  iterates  the  extinction  ratio  for  specific 

%  bit-patterns  and  bit  rates. 


global  B  ns  N  a  I_thresh  I_on  I_off  r  I_p  Power  chrp  p  n  i_p2  t  ... 
init  plot_no; 

%This  routine  iterates  the  value  of  x  that  varies  the  extinction  ratio 


b  =  [.5369]; 
aA  =  [0  1  1  1 
10  0  0 
0  10  1 
10  10]; 

lst_plot  =  0; 

for  j  =  [1:4]; 
a  =  aA(j, 2:4) ; 

init  =  aA(  j , 1) ; 

for  i  =  [1:4] 

B  =  b (1, i) *le9; 

arr  count  =  0; 


for  x  =  [1:. 01:1. 85] 

I  thresh=  33.5e-3; 

I  on  =  1.85*1  thresh;  I _ off  =  x*I _ thresh, 

r= (I_off-I_thresh) / (I_on-I ^thresh) ; 
chrp_mn 

arr_count  =  arr_count  +  1; 
temp  max_chrp  =  max(abs (chrp) ) ; 
temp  r  =  r; 


%Bit  rates  of  interest 
%Bit  patterns  of  interest 


vec (ar recount, 2)  =  temp_max_chrp; 

vec (ar recount , 1 )  —  temper; 


%create  data  arrays  with 
%maximum  excursion  values 


%save  arrays  in  a  workspace 


1st  plot  =  Is t_J?lot  +  1;  %save  arrays  in  a  wori 
savplot  =  [  'save  plot'  num2str (lst_plot )  1  vec  a  B  xmt  ], 
eval  ([savplot]); 
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clear  vec  a  B  init; 


end 

end 
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A-4.  Stage4.m 


Stage4.m  is  a  menu-driven  display  program,  that  allows  the  user  to  specify  specific 
cases  or  a  range  of  cases  to  display. 

Sstage4.m  This  program  allows  the  users  to  call  up  specific  simulation 
%  cases,  and  view  or  print  the  associcated  figure. 

case_in  =  0; 
while  1 

if  case_in  ==  4,  break, end 

clear; 

clg; 

case_in  =  menu (’Which  case  to  Display',  'Specific',  'Range',  ... 

'Print  All',  'Done'); 

SPRINT  ALL  THE  GRAPHICS 
if  case_in  ==  3 
k  =  59; 
while  1 

k=  k+1; 

filek  =  ['case'  int2str(k)]; 
filename  =  [filek  '.mat']; 

if  -exist (filename) , break, end 
eval([*load  '  filename]) 
pause  off 

chrpdsp (t, Power, chrp, p, n, i_p2 , B, a, I_thresh,  I_on,  I_off,  r,  init,  case^no) ; 

print  all  -append  -dps 
pause 
clg 
end 
end 

SPRINT  A  SET  OF  GRAPHIC  IN  A  SPECIFIED  RANGE  WITH  THE  OPTION  OF  PRINTING 
if  case_in  ==  2 

start  =  input (' Starting  Case  #'); 
stop  =  input ('End  Case  #'); 

for  k  =  (start : stop) ; 
filek  =  ['case*  int2str(k)]; 
filename  =  [filek  '.mat']; 
eval(['load  '  filename]) 

chrpdsp (t, Power, chrp,p, n, i_p2, B, a, I _thresh, I_on,  I_of f ,  r,  init,  case_no) ; 
prnt  =  menu ('Print  Out',  'Yes',  'No'); 
orient  tall; 

if  prnt  ==  1, print  -dps, end; 
pause 
end 
end 
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% DISPLAY  A  SPECIFIC  CASE  WITH  THE  OPTION  OF  PRINTING 


if  case_in  —  1 
k  =  input ( ’Case  # ' ) ; 
filek  =  [’case*  int2str(k)]; 
filename  =  [filek  ’.mat'] 
eval ( [ ' load  ’  filename] ) 

chrpdsp (t. Power, chrp,p,n, i_p2, B, a, I_thresh, I_on, I_of f , r, init, case_no) ; 
prnt  =  menu ('Print  Out',  'Yes',  ’No’); 
orient  tall; 

if  prnt  —  1, print  -dps, end; 
pause 
end 

end 
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A-5.  Chrp  dsp.m 


This  m-file  function  plots  the  data  produced  in  chirp. m. 

%chrp_dsp.m  This  function  produces  the  graphic  of  a  specified  case, 

%  ploting  the  data  from  simulation  runs. 

function  [  ]  =chrpdsp (t, Power, chrp,p,n,i_p2,B, a,I_thresh,  ... 
I_on,I_off,r,init, case_no) ; 

orient  tall 
subplot (211) 

y_max  =  max (Power) *le3; 
y_min  =  min (Power) *le3; 

if  max (chrp) *le-9  >  y_max, y_max=max {chrp* le-9) ;, end; 
if  min (chrp) *le-9  <  y_min, y__min=min (chrp*le-9) ; , end; 

y_max  =  round  (yjnax)  +1; 
y_min  =  round ( y_min ) -1 ; 

b=B*le-9; 

plotyy (le9*t, le3*Power, le9*t, le-9*chrp, 1 Frequecy  Chirp  (GHz) 1 ) ; 
xlabel ( 1  time  (ns ) ' ) ;ylabel ( '  Power  (mW) * ) ; 

title ([' Pattern  =>  *num2str (init ) *  1  num2str(a)  '  Bit  rate  =>  1  ... 
num2str(b)  1  Gb/s1]); 

axis([0  max(le9*t)  y_min  y_max] ) ; grid; 
subplot (427) 

plot (le9*t,p, 'b' ) ;  ylabel {' Photon  density,  p');  xlabel ('time  (ns)'); 
title('p(0)  =  1 . 4217el5 ' ) ; axis ( [ 0  le9*max(t)  0  3el5]);grid; 

subplot (428 ) 

plot (le9*t,n, 'b* ) ;  ylabel ( 1  Carrier  density,  n');  xlabel ('time  (ns)'); 
title('n(0)  =  1 . 3978el8 ' ) ; axis ( [ 0  max(le9*t)  1.3el8  1 . 5el8 ] ) ; grid; 

subplot (425) 

plot (le9*t, Ie3*i_p2, 'b' ) ;ylabel ( 'Drive  current  (mA) ' ) ;  xlabel ('time  ... 
(ns)1); 

axis([0  max(le9*t)  0  65]);grid; 

subplot (426) 
axis ( 1  of f '  ) 

text ( 0, . 45, [ ' I_thresh  =  '  num2str  ( I__thresh)  ]  )  ; 
text (0, . 3,  t ' I_of f  =  '  num2str (I_of f ) ] ) ; 
text (0, . 15, [ 'I_on  =  '  num2str(I_on)]); 
text (0, . 9, [ ' r  =  '  num2str (r ) ] ) ; 
text (0,0, [ 'Case  '  num2str (case_no) ] ) ; 

text (0, .75, [ 'Max  chirp  freq  excursion  '  num2str ( le-9*max ( chrp) )  'GHz']); 

text (0, . 6, [ 'Min  chirp  freq  excursion  '  num2str ( le-9*min (chrp) )  '  GHz']); 
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A-6.  Pit  vec.m 


This  program  extracts  data  from  the  data  arrays  produced  in  the  second  simulation 
set.  First,  it  loads  all  the  data  arrays  into  a  matrix  and,  then,  plots  the  extinction  ratio 
versus  maximum  frequency  excursion  curves. 


%plt_vec.m 

% 

% 

% 

% 


This  program  extracts  data  from  the  data  obtained 
in  the  second  simulation  set.  First,  it  loads  all 
the  data  arrays  into  a  matrix  and,  then,  plots  the 
extinction  ratio  versus  maximum  frequency 
excursion  curves. 


clear 

!rm  ext_gr.ps 


for  i  =  [0:3] ; 
for  j  =  [1:4]  ; 

eval  (  ['load  plot'  num2str (i*4+j ) ]  ); 
pltmatx(:,j)  =  vec(:,l); 
pltmaty ( : , j )  =  vec ( : , 2 ) *le-9; 
end 

plot (pltmatx ( : , 1) , pltmaty  ( : , 1) ,pltmatx ( : ,2) , pltmaty ( : ,2) ,  ... 

pltmatx ( : , 3) , pltmaty ( : , 3 ) , pltmatx ( : , 4 ) , pltmaty ( : , 4 ) ) ; 
xlabel (’ Extinction  Ratio'); 
ylabel ( ’Maximum  Frquency  Chirp  (GHz) ' ) ; 

gtext ('0.5  Gb/s ' ) ; 
gtext ( ' 3  Gb/s  * ) ; 
gtext ( ' 6  Gb/s ' ) ; 
gtext ( ' 9  Gb/ s'); 

title ([ 'Maximum  Frequency  Chirp  vs  Extinction  Ratio  (bit  pattern  *  ... 
num2str (init )  '  '  num2str(a)  ')']); 

print  ext_gr  -append 
clg 

end 
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A-7.  Plotyy.m 


This  m-file  is  a  plotting  function  that  allows  both  vertical  sides  of  a  plot  to  have 
separate  labels  and  scales.  This  function  was  provided  from  Mathworks  technical  support 
division. 

function  [ ax, pi , ylab] =plotyy (xl , yl , x2 , y2 , y21abel , options ) 

%  PLOTYY  Plot  graphs  with  Y  tick  labels  on  left  and  right  side. 

%  PLOTYY (xl,yl,x2,y2)  plots  yl  vs.  xl  with  y-axis  labeling 

%  on  left,  and  plots  y2  vs.  x2  with  y-axis  labeling  on  right. 

% 

%  PLOTYY (xl,yl,x2,y2,y21abel)  places  an  optional  y-label  on 

%  the  right  side  of  the  plot,  since  MATLAB  doesn’t  normally 

%  support  right-hand  side  y-labels. 

% 

%  [ax, pi]  =  PLOTYY (xl,yl,x2,y2)  returns  the  handles  of  the  two 

%  axes  created  by  plotyy  in  ax  and  the  handles  of  the  lines 

%  created  by  plotyy  in  pi. 

% 

%  [ax, pi, ylab]  =  PLOTYY (xl, yl, x2, y2, y21abel)  also  returns  the 

%  handle  of  the  text  object  used  for  y21abel  in  ylab. 

% 

%  PLOTYY (xl, yl, x2, y2 , y2 label , options )  allows  for  an  options 

%  vector  as  follows: 

%%  OPTIONS (1)  :  0  for  linear  x-axis,  1  for  log  x-axis 

%  OPTIONS (2)  :  0  for  linear  yl-axis,  1  for  log  yl-axis 

%  OPTIONS (3)  :  0  for  linear  y2-axis,  1  for  log  y2-axis 

%  OPTIONS (4: 5)  :  x-axis  limits,  in  form  [xmin  xmax] 

%  OPTIONS (6: 7)  :  yl-axis  limits,  in  form  [ylmin  ylmax] 

%  OPTIONS (8: 9)  :  y2-axis  limits,  in  form  [y2min  y2max] 

% 

%  Specify  NaN  for  any  option  to  use  default  value. 

% 

%  For  example, 

% 

%  plotyy (xl, yl, x2, y2, ’’,[ 0  1  0  1.5  3  nan  nan  0  10]) 

% 

%  plots  yl  vs.  xl,  y2  vs.  x2,  with  log  yl-axis  scaling  an  x-axis 

%  range  of  1.5  to  3,  the  default  yl-axis  range,  and  a  y2-axis 

%  range  of  0  to  10. 


%Known  bugs : 

%  Due  to  a  bug  in  4.0a  on  Sun  and  SGI  (and  perhaps  others) 

%  setting  the  ylabels  units  to  normalized  (so  they  stick  with  the 
%  axis  as  the  figure  is  resized)  does  not  work  correctly.  In 
%  this  case,  the  ylabel  is  still  drawn,  but  it  won't  follow  the  axis 
%  around  on  the  screen  if  the  figure  is  resized. 


%  constants: 

NUMOPTS  =9;  %  Maximum  number  of  arguments  in  OPTIONS  vector 

01  XLOG  =1;  %  Index  of  Options  vector  for  x-axis  scaling 

OI~ Y1LOG  =2;  %  Index  of  Options  vector  for  yl-axis  scaling 
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01  Y2L0G  =  3;  %  Index  of  Options  vector  for  y2-axis  scaling 
0i_ XMIN  =  4;  %  Index  of  Options  vector  for  x-axis  min  limit 
0I*~XMAX  =5;  %  Index  of  Options  vector  for  x-axis  max  limit 

0I_ Y1MIN  =  6;  %  Index  of  Options  vector  for  yl-axis  min  limit 

0I~Y1MAX  =  7;  %  Index  of  Options  vector  for  yl-axis  max 

0i~ Y2MIN  =8;  %  Index  of  Options  vector  for  y2-axis  min  linut 

0i_ Y2MAX  =9;  %  Index  of  Options  vector  for  y2-axis  max  limit 


%  error  checking: 

if  nargin<4,  , .  . 

error (' PLOTYY  requires  at  least  4  input  arguments.  ), 

elseif  nargin>6, 

error (' PLOTYY  accepts  no  more  than  6  input  arguments.  ), 
elseif  nargout>3, 

error ( 'PLOTYY  returns  a  maximum  of  3  output  arguments.  ), 
end 


%  check  for  options  vector: 

1optionM length ( options ) +1 : NUMOPTS ) =nan*  ones (NUMOPTS-length ( options ) , 1 ) ; 

eoptions=nan*ones(NUMOPTS,l);  %  use  all  default  values 
end 


%  first  calculate  overall  x-axis  range: 
xlmin=min (xl (~isnan (xl ) ) ) ; 
x2min=min(x2(~isnan(x2) ) ) ; 
xlmax=max (xl (-isnan(xl)  )  )  ; 
x2max=max (x2 ( -isnan (x2 ) ) )  ; 

xrange= [min (xlmin, x2min)  max (xlmax, x2max) ] f 


if  ~isnan (options (OI^XMIN) ) , 
xrange= [options (OIJCMIN)  xrange (2) ] ; 
end 

if  -isnan (options (0I_XM AX) ) , 
xrange= [xrange (1)  options (0I_XMAX) ] ; 
end 


%  Set  up  (x2,y2)  plot 

axes2=gca; 

p2=plot (x2 , y2 ) ; 

set (gca, •YDir' , ’rev1 ) ; 

view ( 360~200*eps , -90 )  ; 


first: 


%  This  line  and  the  next  are  the  guts  of  plotyy. 
%  They  allow  the  function  to  display  y  tick 
%  labels  on  the  right  hand  side  of  the  plot. 


%  set  up  (xl,yl)  plot  now: 

%  create  a  new  axis  in  same  location  as  the  old  axis 
axesl=axes (• Position' , get (axes2, ' Pos ' ) ) ; 
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pl=plot  (xl, yl) ; 


%  set  both  axes  to  have  same  x  range: 
set (axesl, 'XLim' , xrange) ; 
set(axes2, 'XLim' , xrange) ; 

%  handle  OPTIONS  vector: 
if  options (OI_XLOG) ==0, 
set (axesl, 'XScale' , ’linear' )  ; 
set (axes2, 'XScale' , 'linear' )  ; 
elseif  options (OI_XLOG) ==1, 
set (axesl, 'XScale' , 'log' ) ; 
set (axes2, 'XScale' , 'log' )  ; 
end 


if  options (0I_Y1L0G) ==0, 
set (axesl, ' YScale’ , 'linear* )  ; 
elseif  options (0I_Y1L0G) ==1, 
set (axesl, ' YScale ' , ' log ' ) ; 
end 


if  options (OI_Y2LOG) ==0, 
set (axes2 , 'YScale' , 'linear' ) ; 
elseif  options (OI_Y2LOG) ==1, 
set (axes 2, ' YScale ' , ' log ' ) ; 
end 


if  -isnan  (options  (0I_Y1MIN)  )  ; 

set (axesl, ' YLim' , get (axesl, ' YLim' ) . * [ 0  1] + [options (0I_Y1MIN)  0] ) 
end 


if  -isnan (options (0I_Y1MAX) ) ; 

set (axesl, ’ YLim' , get (axesl, ' YLim* ) . * [ 1  0] + [ 0  options (0I_Y1MAX) ] ) 
end 


if  -isnan (options (OI_Y2MIN) ) ; 

set (axes 2, 1 YLim* , get (axes 2, *  YLim' )  . * [ 0  1]  +  [options (OI_Y2MIN)  0] ) 
end 


if  -isnan (options (OI_Y2MAX) ) ; 

set (axes2, 'YLim' , get (axes2, 'YLim' ).* [1  0]+[0  options (OI_Y2MAX) ] ) 
end 


%  now  tidy  the  plot  up  a  bit: 

ticklen=get (axesl, 'TickLength' ) ; 

set (axes2, 'TickLength' , [ticklen (1)  ticklen (1) ] ) ; 

set (axes2 , 'TickDir' , get (axesl, 'TickDir ' ) )  ; 


51 


set (axes2, 'XTickLabels ' , ' ' ) ; 
set (axesl, 'Box' , 'off 1 ) ; 


%  If  both  yl  and  y2  are  vectors  (not  matrices)  make 
%  the  color  of  y2  be  the  second  color  in  the  ColorOrder. 
%  (This  is  to  make  the  plot  aesthetic.) 

[ml,nl]“size (yl) ; 

[m2, n2]=size (y2) ; 


if  (ml==l  |  nl==l)  &  (m2==l  |  n2==l), 

ColorOrd=get (axesl, 'ColorOrder1 )  ; 

if  size (ColorOrd, 1) --1,  %  in  case  ColorOrder  only  has  1  color  in  it 
set (p2, 'Color ', ColorOrd (2, : ) ) ; 
end 
end 


%  create  y21abel  if  required: 
if  nargin  >  4, 
if  -is empty (y2 label)  , 
ylab=text (0, 0, y21abel) ; 

save_axesl_units=get (axesl, 'Units' ) ;  %  save  current  state 

set (ylab, 'Units ', 'pixels ') ;  %  convert  everything  to 

set (axesl, 'Units', 'pixels' ) ;  %  same  units 

axeslpos=get (axesl, 'Position' ) ; 

%  get  "typical"  ylabel  offset  from  axes,  add  it  to  width  of  axesl 
%  to  get  y21abel's  x-coordinate : 

typical_ylabel_handle=get (axesl, 'YLabel* ) ; 
typ_ylabel_pos=get ( typical_ylabel  Jiandle, 'Position') ; 
ylabof  f  set=-l*typ__ylabel_jpos  (1)  ; 
ylabxcoordinate-ylabof f set+axeslpos (3) ; 

%  y21abel's  y-coordinate  is  same  as  "typical"  value: 
ylabycoordinate=typ_ylabel_j?os  (2)  ; 

set (ylab, ' Position' , [ylabxcoordinate  ylabycoordinate] ) ; 
set (ylab, ' HorizontalAlignment * , ' center ' ) ; 
set (ylab, ' Rotation' , 90) ; 

%  Due  to  a  bug  in  4.0a  on  Sun  and  SGI  (and  perhaps  others) 

%  setting  the  ylabels  units  to  normalized  (so  they  stick  with  the 
%  axis  as  the  figure  is  resized)  does  not  work  correctly.  In 
%  this  case,  draw  the  ylabel,  but  set  it  so  it  won't  stick  to  axis 
%  position: 

if  -strcmp (version,  ' 4 . 0a ’)  , 

set (ylab, 'Units ' , 'normalized' ) ; 
end 


%  restore  old  'units'  settings: 

set (axesl, 'Units  * , save_axesl_units )  ; 
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%  y2  label  is  usually  clipped  if  using  default  axes  position. 
%  Therefore  shrink  axes  width  by  5% 

%  to  leave  room  for  it: 

axespos=get (axesl, ' Position' ) ; 

set( [axesl,axes2] , ' Position axespos- 10  0  .05  0]); 
end  %  -isempty (y21abel) 
end  %  creation  of  y21abel 


%  return  handles  to  axes  and  lines: 
if  nargout>0, ax=[ axesl ;axes2] ;end 
if  nargout>l,pl= [pl;p2] ;end 
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APPENDIX  B.  BIT-  PATTERN  PLOTS 


This  appendix  contains  the  sixty-four  graphic  produced  in  the  first  group  of  simulations. 
The  sixteen  bit-patterns  are  iterated  against  four  different  bit  rates.  The  patterns  for  0.5  Gb/s  are 
found  on  pages  56-62.  Likewise,  the  patterns  for  3  Gb/s  are  found  on  pages  63-79;  patterns  for  6 
Gb/s  are  found  on  pages  80-86,  and  patterns  for  9  Gb/s  are  found  on  pages  87-103. 
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APPENDIX  C.  EXTINCTION  RATIO  PLOTS 


This  appendix  contains  the  maximum  frequency  excursion  versus  extinction  ratio  plots 
produced  from  the  second  set  of  simulations. 
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Maximum  Frequency  Chirp  vs  Extinction  Ratio  (bit  pattern  1  010) 
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